
% condition: -2 -1 1 2, neg: dark, 2: running
cond=[1 2];

% animals 0: all 1:select
animals=[1:8];

% use active cells only
use_act_only=0;

% normalize to first 2 datapoints
norm_to_baseline=1;

% count active cells instead of taking mean activity
count_act_cells=0;

% %integ is 1 for integration method, 0 for cell counting method
integ=1;
% 
% %threshold 
thresh=1.25;

% using raw data
raw_data_flag=1;

layer=2;

cells_good=[];
cells_total=[];


cnt=0;
cnt_e=[];
cnt_c=[];
act_lp_corr=[];
e_act={[] [] [] [] [] []};
c_act={[] [] [] [] [] []};

tpoints=[-48 -24 12 24 48 72];

e_mat=cell(size(hom_data,2),6);
c_mat=cell(size(hom_data,2),6);

vars_e={'lp_e','mean_templ_e','stim_dur_e','running_act_e'};
vars_c={'lp_c','mean_templ_c','stim_dur_c','running_act_c'};
for ind=1:length(vars_e)
    eval([vars_e{ind} '=-1*ones(length(hom_data),7);']);
end
for ind=1:length(vars_c)
    eval([vars_c{ind} '=-1*ones(length(hom_data),7);']);
end

all_act=[];
for ind=1:length(hom_data)
    act_cell_ind=ones(1,length(sum(hom_data(ind).raw_data(1).sig_act)));
    for knd=1:length(hom_data(ind).raw_data)
        act_cell_ind=act_cell_ind&sum(hom_data(ind).raw_data(knd).sig_act)>10;
    end
    
    for knd=1:length(hom_data(ind).raw_data)
        l_ind=zeros(1,length(hom_data(ind).raw_data(knd).act_code));
        for jnd=1:length(cond)
            try
                l_ind=l_ind+double(hom_data(ind).raw_data(knd).act_code==cond(jnd));
            catch
                l_ind=l_ind+double(hom_data(ind).raw_data(knd).act_code==cond(jnd))';
            end
        end
        l_ind=logical(l_ind);
        l_ind=l_ind(2:4:end);
        cell_ind=ones(1,size(hom_data(ind).raw_data(knd).sig_act,2));

        cell_ind=logical(cell_ind);
        tmp=hom_data(ind).raw_data(knd).sig_act(l_ind,cell_ind);
        tmp=tmp./repmat(median(tmp),size(tmp,1),1);
    end
    
    
    for knd=1:length(hom_data(ind).raw_data)
        curr_ind=find(hom_data(ind).raw_data(knd).timepoint==tpoints);
        l_ind=zeros(1,length(hom_data(ind).raw_data(knd).act_code));
        for jnd=1:length(cond)
            try
                l_ind=l_ind+double(hom_data(ind).raw_data(knd).act_code==cond(jnd));
            catch
                l_ind=l_ind+double(hom_data(ind).raw_data(knd).act_code==cond(jnd))';
            end
        end
        l_ind=logical(l_ind);
        l_ind=l_ind(2:4:end);
        cell_ind=ones(1,size(hom_data(ind).raw_data(knd).sig_act,2));
        cell_ind=logical(cell_ind);
        if use_act_only
            cell_ind=cell_ind&act_cell_ind;
        end
        
        if count_act_cells
            tmp=sum(hom_data(ind).raw_data(knd).sig_act(l_ind,cell_ind))>10/sum(l_ind);
        elseif raw_data_flag
            tmp=hom_data(ind).raw_data(knd).sig_act(l_ind,cell_ind);
            tmp=tmp./repmat(median(tmp),size(tmp,1),1);
            
            if integ==1  %integral method then integ=1, cell sum method integ=0 
                tmp(tmp<thresh)=0;
                tmp=sum(tmp)/sum(l_ind);
            elseif integ==2
                tmp=max(tmp);
            else
                tmp=sum(tmp>thresh)/sum(l_ind);
                all_act(end+1:end+length(tmp))=tmp;
                tmp=tmp>0;
            end
        else
            tmp=sum(hom_data(ind).raw_data(knd).sig_act(l_ind,cell_ind))/sum(l_ind);
        end
        if sum(hom_data(ind).animal_ID==animals)
            if sum(hom_data(ind).layer==layer)
                if hom_data(ind).animal_ID<=4
                    sum(tmp)/length(tmp);
                    cnt_e=[cnt_e ind];
                    e_act{curr_ind}(end+1:end+length(tmp))=tmp;
                    e_mat{ind,curr_ind}(end+1:end+length(tmp))=tmp;
                    mean_templ_e(ind,curr_ind)=mean(hom_data(ind).raw_data(knd).template(:));
                    stim_dur_e(ind,curr_ind)=sum(l_ind);
                    running_act_e(ind,curr_ind)=sum((hom_data(ind).raw_data(knd).act_code)==2)/sum((hom_data(ind).raw_data(knd).act_code)>0);
                else
                    cnt_c=[cnt_c ind];
                    c_act{curr_ind}(end+1:end+length(tmp))=tmp;
                    c_mat{ind,curr_ind}(end+1:end+length(tmp))=tmp;
                    mean_templ_c(ind,curr_ind)=mean(hom_data(ind).raw_data(knd).template(:));
                    stim_dur_c(ind,curr_ind)=sum(l_ind);
                    running_act_c(ind,curr_ind)=sum(abs(hom_data(ind).raw_data(knd).act_code)==2)/sum(abs(hom_data(ind).raw_data(knd).act_code)>0);
                end
                cnt=cnt+1;
            end
        end
    end
end

l_c_mean=[];
l_e_mean=[];
l_e_ste=[];
l_c_ste=[];

for ind=1:length(e_act)
    if norm_to_baseline
        ee=e_act{ind}/nanmean([nanmean(e_act{1}) nanmean(e_act{2})]);
        cc=c_act{ind}/nanmean([nanmean(c_act{1}) nanmean(c_act{2})]);
    else
        ee=e_act{ind};
        cc=c_act{ind};
    end
    
    l_e_mean(ind)=nanmean(ee);
    l_c_mean(ind)=nanmean(cc);
    l_e_ste(ind)=nanstd(ee)/sqrt(sum(~isnan(ee)));
    l_c_ste(ind)=nanstd(cc)/sqrt(sum(~isnan(cc)));
    try
        [~,sig_diff_ttest(ind)]=ttest2(ee,cc);
        sig_diff_rsum(ind)=ranksum(ee,cc);
    catch
    end
end

figure;
clf
hold on
plot(tpoints(1:6),l_e_mean(1:6),'ok','linewidth',2)
plot(tpoints(1:2),l_e_mean(1:2),'k','linewidth',2)
plot(tpoints(3:6),l_e_mean(3:6),'k','linewidth',2)
for ind=1:6
    plot([1 1]*tpoints(ind),[-1 1]*l_e_ste(ind)+l_e_mean(ind),'k')
end

plot(tpoints(1:6),l_c_mean(1:6),'or','linewidth',2)
plot(tpoints(1:2),l_c_mean(1:2),'r','linewidth',2)
plot(tpoints(3:6),l_c_mean(3:6),'r','linewidth',2)
for ind=1:6
    plot([1 1]*tpoints(ind),[-1 1]*l_c_ste(ind)+l_c_mean(ind),'r')
end

set(gca,'xtick',tpoints(1:6))
box off

title_str=({sprintf('Cond: %s - animals; %s',num2str(cond),num2str(animals)) ...
sprintf('ne: %d, nc %d, integ: %d, thresh: %d',length(e_act{3}),length(c_act{3}),integ,thresh)});
title(title_str,'fontweight','bold');
xlabel('Mean over cells')



for ind=1:32
    for knd=1:6
        e_sites(ind,knd)=mean(e_mat{ind,knd});
        c_sites(ind,knd)=mean(c_mat{ind+32,knd});
    end
end
e_sites=e_sites/mean(mean(e_sites(:,1:2)));
m_e_sites=mean(e_sites);
s_e_sites=std(e_sites)/sqrt(size(e_sites,1));
c_sites=c_sites/mean(mean(c_sites(:,1:2)));
m_c_sites=mean(c_sites);
s_c_sites=std(c_sites)/sqrt(size(c_sites,1));


figure;
clf
hold on
plot(tpoints(1:6),m_e_sites(1:6),'ok','linewidth',2)
plot(tpoints(1:2),m_e_sites(1:2),'k','linewidth',2)
plot(tpoints(3:6),m_e_sites(3:6),'k','linewidth',2)
for ind=1:6
    plot([1 1]*tpoints(ind),[-1 1]*s_e_sites(ind)+m_e_sites(ind),'k')
end

plot(tpoints(1:6),m_c_sites(1:6),'or','linewidth',2)
plot(tpoints(1:2),m_c_sites(1:2),'r','linewidth',2)
plot(tpoints(3:6),m_c_sites(3:6),'r','linewidth',2)
for ind=1:6
    plot([1 1]*tpoints(ind),[-1 1]*s_c_sites(ind)+m_c_sites(ind),'r')
end
title(title_str,'fontweight','bold');
xlabel('Mean over regions')

